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Abstract 


Turbomachines for propulsion applications 
operate with many different working fluids and 
flow conditions. The flow may be 
incompressible, such as in the liquid hydrogen 
pump in a rocket engine, or supersonic, such as 
in the turbine which may drive the hydrogen 
pump. Separate codes have traditionally been 
used for incompressible and compressible flow 
solvers. The General Equation Set (GES) 
method can be used to solve both incompressible 
and compressible flows, and it is not restricted to 
perfect gases, as are many compressible-flow 
turbomachinery solvers. An unsteady GES 
turbomachinery flow solver has been developed 
and applied to both air and water flows through 
turbines. It has been shown to be an excellent 
alternative to maintaining two separate codes. 


Roman Symbols: 

e Total internal energy per unit mass 

E, F, G Flux vectors 

J Jacobian of coordinate 

transformation 

k Heat conduction coefficient 

P Pressure 

Q Conserved variables 

Q p Primitive variables 

Re Reynolds number 

T temperature 

u, v, w Physical velocity components 

U, V, W Contravariant velocity components 


Superscripts: 


Pseudo-time step index 
Physical time step index 
viscous 


Nomenclature 


Subscripts: 


Greek Symbols: 


7 . C 

S 

M 

A 

T 

* 

r 


transformed coordinate directions 
Kroneker delta 
viscosity 

turbulent viscosity 

transformed time or shear stress 
(depending on context) 
transformed pseudo-time 


t 

x, y, z 

S-ri-C 


turbulent 

Partial differentiation in physical 
coordinates 

Partial differentiation in 
transformed coordinates 


Other symbols: 


( ) 


Transformed coordinates 
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Introduction 


Turbomachines for propulsion applications 
operate on many different fluids and under a 
wide range of flow conditions. The working 
fluid may be air, liquid or gaseous hydrogen, or 
liquid or gaseous oxygen. The flow may be 
incompressible, such as in the fuel pump in a 
liquid-fuel rocket engine, or supersonic, such as 
in the turbine which may drive the fuel pump. 

Both compressible and incompressible flows are 
governed by the Navier-Stokes equations. 
However, in a nearly-incompressible flow there 
is a great disparity in wave speeds, since the 
speed of sound approaches infinity for a truly 
incompressible fluid. A compressible flow 
solver will encounter numerical stiffness if 
applied to a nearly-incompressible flow, and the 
algorithm will fail. Because of this, it is 
common practice to use one algorithm for 
incompressible flows and a different algorithm 
for compressible flows. 

The General Equation Set (GES) method 1 has 
been developed to handle both compressible and 
incompressible flows. It can be used to solve the 
full, unsteady, 3-D Navier-Stokes equations, and 
with the introduction of a single input flag can 
reduce to the pseudo-compressibility method 
commonly used to solve incompressible flows. 
When used in conjunction with a dual time step, 
it can be used for time accurate simualtions. 

When the Navier-Stokes equations are cast in 
their general form, without applying the perfect 
gas relation, they can be solved for any working 
fluid for which properties are available. The 
GES technique in conjunction with the equations 
cast in this form results in a general flow solver 
applicable to most conditions encountered in 
propulsion turbomachinery applications. 

A code has been written to solve 3-D unsteady 
turbomachinery flows using the GES method. 
Test cases have been run on turbines operating in 
low-speed air (inlet Mach number of 0.07), 
supersonic-inlet air, and water. Results have 
been compared with test data where available, 
and with the results of a well-developed ideal- 
gas, compressible turbomachinery flow 
solver 2 . [Abstract note: The cases run so far are 2- 
D. The final paper will contain 3-D examples.] 


Equations 

The 3-D, unsteady, Navier-Stokes equations may 
be expressed in generalized, curvilinear 
coordinates in dual-time-step form as 

Q, + Q t +E^+P r] +G^ - Re -1 (/s'! + P v + 
( 1 ) 

where 

Q = [ p pu pv pw E f , 

the standard set of conserved variables. The 
fluxes also have their standard definitions. The 

variable T represents the transformed physical 
* 

time step, and T represents the transformed 
pseudo-time step for subiterations. 

Defining the vector of primitive 

variables as 

Q P =[p u v w T ] T , 

the Navier-Stokes equations can be rewritten as 

r,(e,) r +r,(e,),+£ { +£,+G f 2 
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(3) 

In order to pre-condition the equations, p p and 
p T in Eqn. (3) are re-defined as the pseudo- 
properties p' p and p' T , where the p subscript is 
used to indicate preconditioned matrices. The 
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values of the pseudo-properties can be chosen to 
optimize convergence for different types of 
flows, such as incompressible, subsonic 
compressible, supersonic, etc. The 

preconditioned matrix is given by 
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( 4 ) 

Numerical Method 


Eqn (2) will be solved using approximate 
factorization. In the dual-time-step method, the 
solution is iterated to convergence in pseudo- 
time at each physical time step. The physical 
time differencing will have 0(2) accuracy, and 
the pseudo-time differencing will be 0(1). 
(Second-order accuracy is not required for 
pseudo-time since it is driven to convergence 
each step.) 

Since two different time coordinates are being 
used, physical and pseudo, the following 
notation will be used to differentiate between 
numerical time differences: 

A n q) = q> n+lk+ ' -qf 
A k n <p = q> n+l ' k -qf 


and 
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Spatial differencing will be notated by the 
symbol d , with a subscript indicating the 
direction of the differencing. The £ -direction 
inviscid flux derivative is approximated by 




and other inviscid and viscous flux differences 
are approximated in a similar fashion. 

Applying Eqns. (5)-(6) to Eqn. (2) yields 




+ +d rJ F’ 1 + d^G n 
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( 7 ) 


We want to write this equation in terms of 
A k Q p , since that is the variable for which we 
will be solving on the left hand side. Adding and 


subtracting — Q n p +l,k within the parenthesis in 
Eqn.(7) and applying Eqn. (6) yields 


A ,. s-n+\,k+\ l.Jfc 

k<P = < P 

with n indicating the physical time and k 
indicating the pseudo-time. 
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The 0(1) pseudo-time difference term in Eqn. 
(2) is approximated by 
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and the 0(2) physical time difference is 
approximated by 


Applying Eqns. (5), (6), and (8) to Eqn.(7) yields 
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All flux Jacobians are at time level n. The last 
two terms on the right hand side can be 
combined if desired: 
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Since we are solving the overall equation for 
values at pseudo-time level k+ 1, Q p +lk remains 
on the left hand side. 


Before performing approximate factorization, it 
is convenient to define 
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The factorized equation is given by 
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Right Hand Side 


The time term on the right hand side is 
straightforward to implement, since it contains 
no differences. Solutions are required at two 
previous time steps in order to get second order 
temporal accuracy. For the viscous terms, 
standard second-order central differences are 
used, as is common practice. Inviscid terms are 
computed using the approximate Riemann solver 
of Roe 3 . The eigenvalues of the system are 
required for Roe’s scheme. Here, the relevant 

eigenvalues are tose of the matrix ^ p A p , which 
are given by 

4,2,3=^ (12) 


and 



A = p T 0 + pp p h,. 

(14) 

A = p T 0 + pp p hj. 

(15) 

0 = 1- ph p 

(16) 

Also, for later use, we define 

A P =p r 0 + pp p h r 

(17) 
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(18) 


A r ^pJ + ppX 


Similar expressions can be written for the other 
two coordinate directions by substituting V or W 
for U and T] or £ for £ in Eqns.(12) and (13). 



p't(i-pK) 

phf 


Since p T has been set equal to p T , we have 


Left Hand Side 

Ideally, fastest convergence is obtained when the 
same method of differencing is used on both the 
left and right hand sides of Eqn. (11). This is 
possible for low-order schemes since the block 
tridiagonal structure of the equations can be 
maintained. Higher order schemes require larger 
difference stencils, and using them on the left 
hand side precludes the use of block tridiagonal 
solvers. Here, Steger- Warming 4 flux-vector 
splitting is used on the left hand side. This does 
not affect the accuracy of the solution, since the 
left hand side is driven to zero during the 
pseudo-time iterations, and it has been shown to 
be a robust technique. 

Preconditioning Parameters 

Equations (4), (15), (17) and (18) contain two 
yet-to-be-defmed parameters, p p and p, . The 

parameter p T has limited effect on convergence 
rate for most fluids, and is simply set equal 
to p T here. 


1 p T (l-ph p ) 

a p pW 


(21) 


Once the artificial sound speed is defined, the 
preconditioning parameter is known. The 
artificial sound speed, and hence the 
preconditioning parameter, would be expected to 
be a function of the characteristics of flow field. 
Venkateswaran and Merkle 1 have developed a 
method for estimating a p based on the 

calculation of inviscid, viscous, pressure- 
gradient, and unsteady velocity scales. 

The inviscid velocity scale is simply the physical 
velocity, 

V inv = \Ju 2 +v 2 + w 2 

The pressure-gradient velocity scale is given by 



( 22 ) 


The square of the speed of sound can be 
expressed as 


The unsteady velocity scale is given by 
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(19) 


where A is defined in Eqn. (14). If Ain Eqn. 

(19) is replaced by A as defined in Eqn. (15), 
an artificial sound speed can be defined. 
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where l x and l y are characteristic unsteady 
length scales in the x and y directions. 
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Substituting Eqn. (15) into Eqn. (20) and solving 
for p p yields 


The viscous velocity scale is a bit more complex. 
It is shown her in 2-D for simplicity. Variables 
are first defined based on the CFL (Courant 
Friedrichs Lewy) number and VNN (Von 
Neumann Number): 
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and 


P = 


CFL 

VNN 




(23) 


Defining U = £ X U + £ y V + £ Z W (contravariant 
velocity without the time metric) and 
V = Tj x li + TJyV + TJ.w , the viscous velocity 
scale is given by 



V'". = max — , 

0(jB-G)+4±(g+?,+g) 

a i a -v) (PhA 

V{a-v) + ^(rj; +ri;+rj:) ^ A ' 

(24) 

The final velocity scale is then given by 

V P = m in [ max ( V mv , V pgr , V um , V visc ) , a\ 

where a is the speed of sound. 

Code Structure 

The general structure code is based on a well- 
established compressible, 3-D, unsteady 
turbomachinery flow solver. 2 It employs a 
system of overset O-grids and H-grids, with the 
values on the O-H boundaries being updated 
each time step by trilinear interpolation from the 
adjacent grid. The Baldwin-Lomax turbulence 
model 5 is used for turbulence closure. 

Results 

2-D 1-1/2 Stage Turbine 

The first test case is the simulation of the flow 
based on experiments on a 1-1/2 stage Low 
Speed Research Rig (LSRR) turbine at United 
Technologies Research Center 6 . 

The grid dimensions in each of the three blade 
rows are 109x45 in the H-grid and 251x41 in the 
O-grid. The grid is shown in Fig. 6. 


Fig. 6 - Grid, viscous 1-1/2 stage LSRR 
turbine 

Pressure contours are shown in Fig. 7. The 
characteristic accelerations and decelerations are 
seen on the appropriate surfaces of the blades. 
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Fig. 7 - Pressure contours, viscous 1-1/2 stage 
LSRR turbine 

Time-averaged surface pressures are shown in 
Figs. 8-10, comparing the computed values with 
experimental data for each blade row. 
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Fig. 8 - Pressure profile, stator 1, 
viscous 1-1/2 stage LSRR turbine 
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Fig. 9 - Pressure profile, rotor, viscous 
1-1/2 stage LSRR turbine 
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Fig. 10 - Pressure profile, stator 2, 
viscous 1-1/2 stage LSRR turbine 


The match is good for the stators, with some 
discrepancy on the suction surface of the rotor. 

The suction-surface mismatch has been observed 
in a number of previous simulations of this 
configuration by other researchers. 

Figure 1 1 shows entropy contours. The boundary 
layer growth and wake shedding can be seen in 
each blade row. The wake from the first stator 
convects into the rotor passage and takes on the 
characteristic hairpin shape as it becomes 
entrained in the rotor blade surface boundary 
layers. Both the stator- 1 wake and the rotor 
wake then interact with the downstream stator, 
creating a complex pattern of high and low 
entropy fluid. 
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Fig. 11 - Entropy contours, viscous 1-1/2 stage 
LSRR turbine 



The next test case is a simulation of the flow 
through a preliminary geometry of the 
supersonic-inflow Turbine Performance 
Optimization (TPO) turbine 7 . This case was 
chosen to test the shock capturing ability of the 
scheme in order to validate the coding of the 
flux-difference splitting. Two vanes and five 
rotor blades were contained in the simulation. In 
the vanes, each H-grid contained 65x51 points, 
and each O-grid contained 171x31 points. The 
rotor blade H-grids each contained 75x51 points, 
and the O-grids each contained 171x21 points. 
The grid is shown in Fig. 12. 
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Fig. 12 - Grid, TPO turbine 


Pressure contours are shown in Fig. 13. Again, 
the general characteristics of the pressure field 
are as expected. 



Fig. 13 - Pressure contours, TPO turbine 


In order to get a better view of the shock waves, 
Mach number contours are shown in Fig. 14. 
The system of shocks interacting with the 
leading edges of the rotor blades are captured 
quite cleanly. This indicates that the flux- 
difference splitting is operating correctly. 



Fig. 14 - Mach number contours, TPO 
turbine 


Conclusions and Recommendations 

The General Equation Set method has been 
implemented along with general fluid properties 
for turbomachinery applications. A series of 
subsonic and supersonic test cases have been 
run, and the code has been demonstrated to 
perform well. 

The dual time step algorithm requires 
specification of a parameter (the ratio of pseudo- 
time-step to physical rime step), and the choice 
will affect the convergence rate. It is 
recommended that a series of simulations be 
performed on typical configurations using a 
range of values for this parameter. This will help 
to deduce a range of the parameter’s value to 
optimize convergence rates. This will be 
particularly for computing incompressible and 
nearly-incompressible flows. 

The code has been constructed in a modular 
fashion to make future enhancements as easy as 
possible. One such enhancement would be to 
develop a series of routines to compute the 
properties for a variety of fluids encountered in 
rocket propulsion turbomachinery (liquid and 
gaseous hydrogen, oxygen, etc.). 

Additional Results to be Shown for 
Final Paper 

The final paper will show results for water flow 
(presently underway), and also at least one 3-D 
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result. If necessary, the text and equations in this 
abstract will be cut down for brevity. 
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